<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<title>AMBHAS: /home/tomer/svn/ambhas/ambhas/krige.py Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<script type="text/javascript">
  $(document).ready(initResizable);
</script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
  $(document).ready(function() { searchBox.OnSelectItem(0); });
</script>

</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">AMBHAS
   
   </div>
   
  </td>
  
  
  
   
   <td>        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.png"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
          </span>
        </div>
</td>
   
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.6.1 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
</div>
<div id="side-nav" class="ui-resizable side-nav-resizable">
  <div id="nav-tree">
    <div id="nav-tree-contents">
    </div>
  </div>
  <div id="splitbar" style="-moz-user-select:none;" 
       class="ui-resizable-handle">
  </div>
</div>
<script type="text/javascript">
  initNavTree('krige_8py.html','');
</script>
<div id="doc-content">
<div class="header">
  <div class="headertitle">
<div class="title">krige.py</div>  </div>
</div><!--header-->
<div class="contents">
<a href="krige_8py.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a><a class="code" href="namespaceambhas_1_1krige.html">00001</a> <span class="comment"># -*- coding: utf-8 -*-</span>
<a name="l00002"></a>00002 <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00003"></a>00003 <span class="stringliteral">Created on Thu Jun  9 17:55:54 2011</span>
<a name="l00004"></a>00004 <span class="stringliteral"></span>
<a name="l00005"></a>00005 <span class="stringliteral">@author: Sat Kumar Tomer</span>
<a name="l00006"></a>00006 <span class="stringliteral">@website: www.ambhas.com</span>
<a name="l00007"></a>00007 <span class="stringliteral">@email: satkumartomer@gmail.com</span>
<a name="l00008"></a>00008 <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00009"></a>00009 
<a name="l00010"></a>00010 <span class="comment"># import required modules</span>
<a name="l00011"></a>00011 <span class="keyword">import</span> numpy <span class="keyword">as</span> np
<a name="l00012"></a>00012 <span class="keyword">import</span> matplotlib.pylab <span class="keyword">as</span> plt
<a name="l00013"></a>00013 
<a name="l00014"></a>00014 
<a name="l00015"></a><a class="code" href="classambhas_1_1krige_1_1OK.html">00015</a> <span class="keyword">class </span><a class="code" href="classambhas_1_1krige_1_1OK.html">OK</a>:
<a name="l00016"></a>00016     <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00017"></a>00017 <span class="stringliteral">    This performs the ordinary kriging</span>
<a name="l00018"></a>00018 <span class="stringliteral">    Input:</span>
<a name="l00019"></a>00019 <span class="stringliteral">        x: x vector of location</span>
<a name="l00020"></a>00020 <span class="stringliteral">        Y: y vector of location</span>
<a name="l00021"></a>00021 <span class="stringliteral">        z: data vector at location (x,y)</span>
<a name="l00022"></a>00022 <span class="stringliteral">    </span>
<a name="l00023"></a>00023 <span class="stringliteral">    Output:</span>
<a name="l00024"></a>00024 <span class="stringliteral">        None</span>
<a name="l00025"></a>00025 <span class="stringliteral">        </span>
<a name="l00026"></a>00026 <span class="stringliteral">    Methods:</span>
<a name="l00027"></a>00027 <span class="stringliteral">        variogram: estimate the variogram</span>
<a name="l00028"></a>00028 <span class="stringliteral">        </span>
<a name="l00029"></a>00029 <span class="stringliteral">    &quot;&quot;&quot;</span>
<a name="l00030"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a12d3a9f59c4b6637aad1dde4e8b7bfdd">00030</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1krige_1_1OK.html#a1a4f47e40d58c335a4d4842a2f996d61">__init__</a>(self,x,y,z):
<a name="l00031"></a>00031         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a> = x.flatten()
<a name="l00032"></a>00032         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#addb0c60b4f851d1667f59fcb0340fe4f">y</a> = y.flatten()
<a name="l00033"></a>00033         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a12d3a9f59c4b6637aad1dde4e8b7bfdd">z</a> = z.flatten()
<a name="l00034"></a>00034     
<a name="l00035"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a70d979f4120eb6f610a171bc0edd0592">00035</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1krige_1_1OK.html#a70d979f4120eb6f610a171bc0edd0592">variogram</a>(self, var_type=&#39;averaged&#39;, n_lag=9):
<a name="l00036"></a>00036         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00037"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a2317fe0fde2866e1e986e2fd802cd6ab">00037</a> <span class="stringliteral">        var_type: averaged or scattered</span>
<a name="l00038"></a>00038 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00039"></a>00039         
<a name="l00040"></a>00040         x = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a>
<a name="l00041"></a>00041         y = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#addb0c60b4f851d1667f59fcb0340fe4f">y</a>
<a name="l00042"></a>00042         z = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a12d3a9f59c4b6637aad1dde4e8b7bfdd">z</a>
<a name="l00043"></a>00043         <span class="comment"># make the meshgrid</span>
<a name="l00044"></a>00044         X1,X2 = np.meshgrid(x,x) 
<a name="l00045"></a>00045         Y1,Y2 = np.meshgrid(y,y)
<a name="l00046"></a>00046         Z1,Z2 = np.meshgrid(z,z)
<a name="l00047"></a>00047         
<a name="l00048"></a>00048         D = np.sqrt((X1 - X2)**2 + (Y1 - Y2)**2)
<a name="l00049"></a>00049         
<a name="l00050"></a>00050         G = 0.5*(Z1 - Z2)**2
<a name="l00051"></a>00051         indx = range(len(z))
<a name="l00052"></a>00052         C,R = np.meshgrid(indx,indx)
<a name="l00053"></a>00053         G = G[R&gt;C]
<a name="l00054"></a>00054         
<a name="l00055"></a>00055         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a2317fe0fde2866e1e986e2fd802cd6ab">D</a> = D
<a name="l00056"></a>00056         DI = D[R &gt; C]
<a name="l00057"></a>00057         
<a name="l00058"></a>00058         <span class="comment"># group the variogram</span>
<a name="l00059"></a>00059         <span class="comment"># the group are formed based on the equal number of bin</span>
<a name="l00060"></a>00060         total_n = len(DI)
<a name="l00061"></a>00061         group_n = int(total_n/n_lag)
<a name="l00062"></a>00062         sor_i = np.argsort(DI)[::-1]
<a name="l00063"></a>00063         
<a name="l00064"></a>00064         DE = np.empty(n_lag)
<a name="l00065"></a>00065         GE = np.empty(n_lag)
<a name="l00066"></a>00066         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(n_lag):
<a name="l00067"></a>00067             <span class="keywordflow">if</span> i&lt;n_lag-1:
<a name="l00068"></a>00068                 DE[i] = DI[sor_i[group_n*i:group_n*(i+1)]].mean()
<a name="l00069"></a>00069                 GE[i] = G[sor_i[group_n*i:group_n*(i+1)]].mean()
<a name="l00070"></a>00070                 
<a name="l00071"></a>00071             <span class="keywordflow">else</span>:
<a name="l00072"></a>00072                 DE[i] = DI[sor_i[group_n*i:]].mean()
<a name="l00073"></a>00073                 GE[i] = G[sor_i[group_n*i:]].mean()
<a name="l00074"></a>00074             
<a name="l00075"></a>00075         <span class="keywordflow">if</span> var_type == <span class="stringliteral">&#39;scattered&#39;</span>:
<a name="l00076"></a>00076             <span class="keywordflow">return</span> DI,G      
<a name="l00077"></a>00077         <span class="keywordflow">elif</span> var_type == <span class="stringliteral">&#39;averaged&#39;</span>:
<a name="l00078"></a>00078             <span class="keywordflow">return</span> DE,GE
<a name="l00079"></a>00079         <span class="keywordflow">else</span>:
<a name="l00080"></a>00080             <span class="keywordflow">raise</span> ValueError(<span class="stringliteral">&#39;var_type should be either averaged or scatter&#39;</span>)
<a name="l00081"></a>00081         
<a name="l00082"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a72309c8a1f65885c87d832805e4ba629">00082</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1krige_1_1OK.html#a72309c8a1f65885c87d832805e4ba629">vario_model</a>(self, lags, model_par, model_type=&#39;linear&#39;):
<a name="l00083"></a>00083         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00084"></a>00084 <span class="stringliteral">        Input:</span>
<a name="l00085"></a>00085 <span class="stringliteral">            model_type : the type of variogram model </span>
<a name="l00086"></a>00086 <span class="stringliteral">                             spherical</span>
<a name="l00087"></a>00087 <span class="stringliteral">                             linear</span>
<a name="l00088"></a>00088 <span class="stringliteral">                             exponential</span>
<a name="l00089"></a>00089 <span class="stringliteral">            model_par:  parameters of variogram model</span>
<a name="l00090"></a>00090 <span class="stringliteral">                        this should be a dictionary </span>
<a name="l00091"></a>00091 <span class="stringliteral">                        e.g. for shperical and exponential</span>
<a name="l00092"></a>00092 <span class="stringliteral">                            model_par = {&#39;nugget&#39;:0, &#39;range&#39;:1, &#39;sill&#39;:1}</span>
<a name="l00093"></a>00093 <span class="stringliteral">                        for linear</span>
<a name="l00094"></a>00094 <span class="stringliteral">                            model_par = {&#39;nugget&#39;:0, &#39;slope&#39;:1}</span>
<a name="l00095"></a>00095 <span class="stringliteral">        Output:</span>
<a name="l00096"></a>00096 <span class="stringliteral">            G:  The fitted variogram model</span>
<a name="l00097"></a>00097 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00098"></a>00098         
<a name="l00099"></a>00099         <span class="keywordflow">if</span> model_type == <span class="stringliteral">&#39;spherical&#39;</span>:
<a name="l00100"></a>00100             n = model_par[<span class="stringliteral">&#39;nugget&#39;</span>]
<a name="l00101"></a>00101             r = model_par[<span class="stringliteral">&#39;range&#39;</span>]
<a name="l00102"></a>00102             s = model_par[<span class="stringliteral">&#39;sill&#39;</span>]
<a name="l00103"></a>00103             l = lags
<a name="l00104"></a>00104             G = n + (s*(1.5*l/r - 0.5*(l/r)**3)*(l&lt;=r) + s*(l&gt;r))
<a name="l00105"></a>00105         
<a name="l00106"></a>00106         <span class="keywordflow">elif</span> model_type == <span class="stringliteral">&#39;linear&#39;</span>:
<a name="l00107"></a>00107             n = model_par[<span class="stringliteral">&#39;nugget&#39;</span>]
<a name="l00108"></a>00108             s = model_par[<span class="stringliteral">&#39;slope&#39;</span>]
<a name="l00109"></a>00109             l = lags
<a name="l00110"></a>00110             G = n + s*l
<a name="l00111"></a>00111         
<a name="l00112"></a>00112         <span class="keywordflow">elif</span> model_type == <span class="stringliteral">&#39;exponential&#39;</span>:
<a name="l00113"></a>00113             n = model_par[<span class="stringliteral">&#39;nugget&#39;</span>]
<a name="l00114"></a>00114             r = model_par[<span class="stringliteral">&#39;range&#39;</span>]
<a name="l00115"></a>00115             s = model_par[<span class="stringliteral">&#39;sill&#39;</span>]
<a name="l00116"></a>00116             l = lags
<a name="l00117"></a>00117             G = n + s*(1 - np.exp(-3*l/r))
<a name="l00118"></a>00118         
<a name="l00119"></a>00119         <span class="keywordflow">else</span>:
<a name="l00120"></a>00120             <span class="keywordflow">raise</span> ValueError(<span class="stringliteral">&#39;model_type should be spherical or linear or exponential&#39;</span>)
<a name="l00121"></a>00121             
<a name="l00122"></a>00122         <span class="keywordflow">return</span> G
<a name="l00123"></a>00123 
<a name="l00124"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a51b9c922a7853872fab765745c9dccdf">00124</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1krige_1_1OK.html#a51b9c922a7853872fab765745c9dccdf">int_vario</a>(self, Xg, Yg, model_par, model_type):
<a name="l00125"></a>00125         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00126"></a>00126 <span class="stringliteral">        this computes the integral of the variogram over a square</span>
<a name="l00127"></a>00127 <span class="stringliteral">        using the Monte Carlo integration method</span>
<a name="l00128"></a>00128 <span class="stringliteral">        </span>
<a name="l00129"></a>00129 <span class="stringliteral">        this works only for two dimensional grid</span>
<a name="l00130"></a>00130 <span class="stringliteral">        </span>
<a name="l00131"></a>00131 <span class="stringliteral">        Input:</span>
<a name="l00132"></a>00132 <span class="stringliteral">            Xg:     x location where krigged data is required</span>
<a name="l00133"></a>00133 <span class="stringliteral">            Yg:     y location whre kirgged data is required</span>
<a name="l00134"></a>00134 <span class="stringliteral">            model_par: see the vario_model</span>
<a name="l00135"></a>00135 <span class="stringliteral">            model_type: see the vario_model</span>
<a name="l00136"></a>00136 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00137"></a>00137         avg_vario = np.empty((len(self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a>), (len(Xg)-1)*(len(Yg)-1)))
<a name="l00138"></a>00138         <span class="keywordflow">for</span> k <span class="keywordflow">in</span> range(len(self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a>)):
<a name="l00139"></a>00139             
<a name="l00140"></a>00140             avg_vario_ens = np.empty((len(Xg)-1, len(Yg)-1))
<a name="l00141"></a>00141             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(len(Xg)-1):
<a name="l00142"></a>00142                 <span class="keywordflow">for</span> j <span class="keywordflow">in</span> range(len(Yg)-1):
<a name="l00143"></a>00143                     Xg_rand = Xg[i]+np.random.rand(10)*(Xg[i+1]-Xg[i])
<a name="l00144"></a>00144                     Yg_rand = Yg[j]+np.random.rand(10)*(Yg[j+1]-Yg[j])    
<a name="l00145"></a>00145 
<a name="l00146"></a>00146                     DOR = ((self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a>[k] - Xg_rand)**2 + (self.<a class="code" href="classambhas_1_1krige_1_1OK.html#addb0c60b4f851d1667f59fcb0340fe4f">y</a>[k] - Yg_rand)**2)**0.5
<a name="l00147"></a>00147                     avg_vario_ens[i,j] = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a72309c8a1f65885c87d832805e4ba629">vario_model</a>(DOR, model_par, model_type).mean()
<a name="l00148"></a>00148             avg_vario[k,:] = avg_vario_ens.flatten()
<a name="l00149"></a>00149         <span class="keywordflow">return</span> avg_vario
<a name="l00150"></a>00150     
<a name="l00151"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a076f9071674fbbea73df1ec0a746e2cc">00151</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1krige_1_1OK.html#a076f9071674fbbea73df1ec0a746e2cc">krige</a>(self, Xg, Yg, model_par, model_type):
<a name="l00152"></a>00152         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00153"></a>00153 <span class="stringliteral">        Input:</span>
<a name="l00154"></a>00154 <span class="stringliteral">            Xg:     x location where krigged data is required</span>
<a name="l00155"></a>00155 <span class="stringliteral">            Yg:     y location whre kirgged data is required</span>
<a name="l00156"></a>00156 <span class="stringliteral">            model_par: see the vario_model</span>
<a name="l00157"></a>00157 <span class="stringliteral">            model_type: see the vario_model</span>
<a name="l00158"></a>00158 <span class="stringliteral">            </span>
<a name="l00159"></a>00159 <span class="stringliteral">        Attributes:</span>
<a name="l00160"></a>00160 <span class="stringliteral">            self.Zg : krigged data</span>
<a name="l00161"></a>00161 <span class="stringliteral">            self.s2_k = variance in the data</span>
<a name="l00162"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a25442a099a818290d039b304561c3f05">00162</a> <span class="stringliteral">                </span>
<a name="l00163"></a>00163 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00164"></a>00164         
<a name="l00165"></a>00165         <span class="comment"># set up the Gmod matrix </span>
<a name="l00166"></a>00166         n = len(self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a>)
<a name="l00167"></a>00167         Gmod = np.empty((n+1,n+1))
<a name="l00168"></a>00168         Gmod[:n, :n] = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a72309c8a1f65885c87d832805e4ba629">vario_model</a>(self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a2317fe0fde2866e1e986e2fd802cd6ab">D</a>, model_par, model_type)
<a name="l00169"></a>00169                 
<a name="l00170"></a>00170         Gmod[:,n] = 1
<a name="l00171"></a>00171         Gmod[n,:] = 1
<a name="l00172"></a>00172         Gmod[n,n] = 0
<a name="l00173"></a>00173 
<a name="l00174"></a>00174         Gmod = np.matrix(Gmod)      
<a name="l00175"></a>00175         
<a name="l00176"></a>00176         <span class="comment"># inverse of Gmod</span>
<a name="l00177"></a>00177         Ginv = Gmod.I
<a name="l00178"></a>00178 
<a name="l00179"></a>00179         Xg = Xg.flatten()
<a name="l00180"></a>00180         Yg = Yg.flatten()        
<a name="l00181"></a>00181         Zg = np.empty(Xg.shape)
<a name="l00182"></a>00182         s2_k = np.empty(Xg.shape)
<a name="l00183"></a>00183         
<a name="l00184"></a>00184         <span class="keywordflow">for</span> k <span class="keywordflow">in</span> range(len(Xg)):
<a name="l00185"></a>00185             
<a name="l00186"></a>00186             DOR = ((self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a> - Xg[k])**2 + (self.<a class="code" href="classambhas_1_1krige_1_1OK.html#addb0c60b4f851d1667f59fcb0340fe4f">y</a> - Yg[k])**2)**0.5
<a name="l00187"></a>00187             GR = np.empty((n+1,1))
<a name="l00188"></a>00188             
<a name="l00189"></a>00189             GR[:n,0] = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a72309c8a1f65885c87d832805e4ba629">vario_model</a>(DOR, model_par, model_type)
<a name="l00190"></a>00190             
<a name="l00191"></a>00191             GR[n,0] = 1
<a name="l00192"></a>00192             E = np.array(Ginv * GR )
<a name="l00193"></a>00193             Zg[k] = np.sum(E[:n,0]*self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a12d3a9f59c4b6637aad1dde4e8b7bfdd">z</a>)
<a name="l00194"></a>00194             s2_k[k] = np.sum(E[:n,0]*GR[:n,0])+ E[n, 0]
<a name="l00195"></a>00195         
<a name="l00196"></a>00196         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a25442a099a818290d039b304561c3f05">Zg</a> = Zg
<a name="l00197"></a>00197         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a6da5cd342069a87c010de2af4eb1ce38">s2_k</a> = s2_k
<a name="l00198"></a>00198         
<a name="l00199"></a><a class="code" href="classambhas_1_1krige_1_1OK.html#a046ca0bd7a04c619ec9f5727d867ecd2">00199</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1krige_1_1OK.html#a046ca0bd7a04c619ec9f5727d867ecd2">block_krige</a>(self, Xg, Yg, model_par, model_type):
<a name="l00200"></a>00200         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00201"></a>00201 <span class="stringliteral">        Input:</span>
<a name="l00202"></a>00202 <span class="stringliteral">            Xg:     x location where krigged data is required</span>
<a name="l00203"></a>00203 <span class="stringliteral">            Yg:     y location whre krigged data is required</span>
<a name="l00204"></a>00204 <span class="stringliteral">            model_par: see the vario_model</span>
<a name="l00205"></a>00205 <span class="stringliteral">            model_type: see the vario_model</span>
<a name="l00206"></a>00206 <span class="stringliteral">            </span>
<a name="l00207"></a>00207 <span class="stringliteral">        Attributes:</span>
<a name="l00208"></a>00208 <span class="stringliteral">            self.Zg : krigged data</span>
<a name="l00209"></a>00209 <span class="stringliteral">            self.s2_k = variance in the data</span>
<a name="l00210"></a>00210 <span class="stringliteral">                </span>
<a name="l00211"></a>00211 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00212"></a>00212         
<a name="l00213"></a>00213         <span class="comment"># set up the Gmod matrix </span>
<a name="l00214"></a>00214         n = len(self.<a class="code" href="classambhas_1_1krige_1_1OK.html#adfe2a29b2df4daac8a278ef6253e0c94">x</a>)
<a name="l00215"></a>00215         Gmod = np.empty((n+1,n+1))
<a name="l00216"></a>00216         Gmod[:n, :n] = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a72309c8a1f65885c87d832805e4ba629">vario_model</a>(self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a2317fe0fde2866e1e986e2fd802cd6ab">D</a>, model_par, model_type)
<a name="l00217"></a>00217                 
<a name="l00218"></a>00218         Gmod[:,n] = 1
<a name="l00219"></a>00219         Gmod[n,:] = 1
<a name="l00220"></a>00220         Gmod[n,n] = 0
<a name="l00221"></a>00221 
<a name="l00222"></a>00222         Gmod = np.matrix(Gmod)      
<a name="l00223"></a>00223         
<a name="l00224"></a>00224         <span class="comment"># inverse of Gmod</span>
<a name="l00225"></a>00225         Ginv = Gmod.I
<a name="l00226"></a>00226 
<a name="l00227"></a>00227         Xg = Xg.flatten()
<a name="l00228"></a>00228         Yg = Yg.flatten()        
<a name="l00229"></a>00229         
<a name="l00230"></a>00230      
<a name="l00231"></a>00231         avg_vario = self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a51b9c922a7853872fab765745c9dccdf">int_vario</a>(Xg, Yg, model_par, model_type)
<a name="l00232"></a>00232         Zg = np.empty(avg_vario.shape[1])
<a name="l00233"></a>00233         s2_k = np.empty(avg_vario.shape[1])
<a name="l00234"></a>00234         
<a name="l00235"></a>00235         <span class="keywordflow">for</span> k <span class="keywordflow">in</span> range(avg_vario.shape[1]):
<a name="l00236"></a>00236             
<a name="l00237"></a>00237             GR = np.empty((n+1,1))
<a name="l00238"></a>00238             GR[:n,0] = avg_vario[:,k]
<a name="l00239"></a>00239             GR[n,0] = 1
<a name="l00240"></a>00240             E = np.array(Ginv * GR )
<a name="l00241"></a>00241             Zg[k] = np.sum(E[:n,0]*self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a12d3a9f59c4b6637aad1dde4e8b7bfdd">z</a>)
<a name="l00242"></a>00242             s2_k[k] = np.sum(E[:n,0]*GR[:n,0])+ E[n, 0]
<a name="l00243"></a>00243         
<a name="l00244"></a>00244         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a25442a099a818290d039b304561c3f05">Zg</a> = Zg.reshape(len(Xg)-1, len(Yg)-1)
<a name="l00245"></a>00245         self.<a class="code" href="classambhas_1_1krige_1_1OK.html#a6da5cd342069a87c010de2af4eb1ce38">s2_k</a> = s2_k.reshape(len(Xg)-1, len(Yg)-1)
<a name="l00246"></a>00246             
<a name="l00247"></a>00247 <span class="keywordflow">if</span> __name__ == <span class="stringliteral">&quot;__main__&quot;</span>:          
<a name="l00248"></a>00248     <span class="comment"># generate some sythetic data</span>
<a name="l00249"></a><a class="code" href="namespaceambhas_1_1krige.html#adc16e92eaa5a77b851e4c33cf80be64f">00249</a>     x = np.random.rand(20)
<a name="l00250"></a><a class="code" href="namespaceambhas_1_1krige.html#ae890db33bed8b2d3e4f2cd67f40b3630">00250</a>     y = np.random.rand(20)
<a name="l00251"></a><a class="code" href="namespaceambhas_1_1krige.html#a05974f3d5547c0bba87c55b759124824">00251</a>     z = 0.0*np.random.normal(size=20)+x+y
<a name="l00252"></a>00252     
<a name="l00253"></a><a class="code" href="namespaceambhas_1_1krige.html#a7cb2cd35072459de5278d74f36718f3a">00253</a>     foo = <a class="code" href="classambhas_1_1krige_1_1OK.html">OK</a>(x,y,z)
<a name="l00254"></a>00254     <span class="comment">#ax,ay = foo.variogram(&#39;scattered&#39;)</span>
<a name="l00255"></a>00255     ax,ay = foo.variogram()
<a name="l00256"></a>00256     
<a name="l00257"></a>00257     plt.plot(ax,ay,<span class="stringliteral">&#39;ro&#39;</span>)
<a name="l00258"></a>00258     
<a name="l00259"></a><a class="code" href="namespaceambhas_1_1krige.html#a5383ff41eb19e27c42dcaa6bc4873dd2">00259</a>     lags = np.linspace(0,5)
<a name="l00260"></a><a class="code" href="namespaceambhas_1_1krige.html#a0b00f5f8782630cad0eebbad3ee4ede3">00260</a>     model_par = {}
<a name="l00261"></a>00261     model_par[<span class="stringliteral">&#39;nugget&#39;</span>] = 0
<a name="l00262"></a>00262     model_par[<span class="stringliteral">&#39;range&#39;</span>] = 1
<a name="l00263"></a>00263     model_par[<span class="stringliteral">&#39;sill&#39;</span>] = 2.0
<a name="l00264"></a>00264     
<a name="l00265"></a><a class="code" href="namespaceambhas_1_1krige.html#a68d92bac7990a6eee5bb9d0b76adddbc">00265</a>     G = foo.vario_model(lags, model_par, model_type = <span class="stringliteral">&#39;exponential&#39;</span>)
<a name="l00266"></a>00266     plt.plot(lags, G, <span class="stringliteral">&#39;k&#39;</span>)
<a name="l00267"></a>00267     plt.show()
<a name="l00268"></a>00268     
<a name="l00269"></a><a class="code" href="namespaceambhas_1_1krige.html#aae8f3009377fa08789753dff30bc5119">00269</a>     Rx = np.linspace(-1,1,1050)
<a name="l00270"></a><a class="code" href="namespaceambhas_1_1krige.html#ac474e70ee53cf6c4a238445daa7418db">00270</a>     Ry = np.linspace(0,1,750)
<a name="l00271"></a>00271     XI,YI = np.meshgrid(Rx,Ry)
<a name="l00272"></a>00272     foo.krige(XI, YI, model_par, <span class="stringliteral">&#39;exponential&#39;</span>)
<a name="l00273"></a>00273     
<a name="l00274"></a>00274     plt.matshow(foo.Zg.reshape(750,1050))
<a name="l00275"></a>00275     plt.show()
<a name="l00276"></a>00276     
<a name="l00277"></a>00277 <span class="comment">#    # block kriging</span>
<a name="l00278"></a>00278 <span class="comment">#    xg = np.linspace(0,1,5)</span>
<a name="l00279"></a>00279 <span class="comment">#    yg = np.linspace(0,1,8)</span>
<a name="l00280"></a>00280 <span class="comment">#    foo.block_krige(xg, yg, model_par, model_type = &#39;exponential&#39;)</span>
<a name="l00281"></a>00281 <span class="comment">#    plt.imshow(foo.s2_k, extent=(0,1,0,1))</span>
<a name="l00282"></a>00282 <span class="comment">#    plt.imshow(foo.Zg, extent=(0,1,0,1))</span>
<a name="l00283"></a>00283 <span class="comment">#    plt.matshow(foo.Zg)</span>
<a name="l00284"></a>00284 <span class="comment">#    plt.matshow(foo.s2_k)</span>
<a name="l00285"></a>00285 <span class="comment">#    plt.colorbar()</span>
<a name="l00286"></a>00286 <span class="comment">#    plt.plot(x,y, &#39;ro&#39;)</span>
<a name="l00287"></a>00287 <span class="comment">#    plt.show()</span>
</pre></div></div><!-- contents -->
</div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Namespaces</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark">&#160;</span>Variables</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

  <div id="nav-path" class="navpath">
    <ul>
      <li class="navelem"><a class="el" href="krige_8py.html">krige.py</a>      </li>

    <li class="footer">Generated on Sat Jul 21 2012 12:26:08 for AMBHAS by
    <a href="http://www.doxygen.org/index.html">
    <img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.6.1 </li>
   </ul>
 </div>


</body>
</html>
